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Abstract 

We propose a simple mathematical model (a coupled system of nonlinear 
ODEs) able to capture dynamical effects produced by adding charcoal to 
fertile soils. Our main aim is to understand to which extent charcoal (in its 
biochar form) is able to lock up carbon in soils. Our results are preliminary 
in the sense that we do not actually solve the CO2 sequestration problem, but 
we do set up a modeling framework in which this can be tackled by means 
of mathematical tools. 

We show that our model is well-posed and has interesting steady states. 
Depending on the reference parameter range and chosen time scale, numerical 
simulations suggest that adding charcoal postpones the release of CO2 for a 
large variety of soils. 

Keywords: Modeling chemical kinetics in fertile soils, Solvability of a 
nonlinear ODE system, Equilibria and steady states, Simulation, Biochar, 
CO2 sequestration 



1. Introduction 

Global warming is an increasingly important issue for mankind. It seems 
that it is no longer enough to reduce CO2 emissions; one also needs to remove 
CO2 from the atmosphere (carbon sequestration). In his Nature paper [I], J. 
Lehmann argues that locking carbon up in soil makes more sense than storing 
it in plants and trees that eventually decompose, but can this idea work on a 
large timescale? A large community of scientists (mostly biologists, chemists, 
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soil engineers) started to support such ideas and tried, with their experi- 
mental means, to explore the sustainability of adding charcoal (biochar) to 
soils; compare for instance [21 El SJ and see also the review paper [5]. For 
more information on this research trend, often called the Biochar project, 
we refer the reader also to the sites www.biochar-international.org and 
http : //en. wikipedia. org/wiki/Biochar. The Biochar project]^ brings clear 
advantages^] (e.g. reduces soil greenhouse gas emissions, improves both water 
holding and nutrient holding capacities, improves environment for soil life, 
doesn't alter the carbon/nitrogen ratio, reduces soil acidity, remove pollu- 
tants), but is it a secure solution? Is it a permanent one? What about the 
possible negative effects like charcoal increases soil fertility and so increases 
the microbe population, which finally can lead to an increase in the rate at 
which the natural soil carbon is broken down and then released as CO2? 
In spite of the intense current experimental research, there is no conclusive 
evidence yet about whether putting charcoal in soil is a good idea, therefore 
our interest. 

In our opinion, Lehmann's question (see loc. cit. in pQ) can be trans- 
lated into mathematical terms as follows: What is the large time behavior 
of the complex dynamical system (including transport, soil geometry and 
chemistry) provoked by adding charcoal? The major issue is the complexity 
of the situation - it is a priori not clear what a good model is to capture 
the effect of charcoal on CO2 emissions. This is the place where we wish to 
contribute. 

In particular, note that charcoal (or biochar) is characterized by a very 
special porous structure (see Figure [I]), which is responsible for the high 
retention of water, dissolved organic nutrients, and even of pollutants such 
as hydrocarbons and pesticides^} On top of this, the chemistry of fertile soils 
is rather complex and precise characterizations of the microbial evolution are 
not available. Furthermore, describing the transport of water together with 
nutrients, phenolics, pollutants (etc.) requires a good understanding of the 
heterogeneities of the soils. 



1 Biochar := The idea of trapping carbon in soil for longer by storing it in the form of 
charcoal. 

2 Note also the additional advantage of producing energy by burning organic matter to 
make charcoal. 

3 As a direct consequence of this fact, there are several situations when the soil fertility 
has increased after charcoal addition. 
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Figure 1: Multiscale geometry of biochar (left: macro, right: micro). This is the place 
where nutrients, phenolics etc. undergo adsorption and desorption. 



Within this framework we treat a spatially homogeneous soil. Hence we 
avoid the aforementioned complications and propose the simplest mathe- 
matical model that is able to capture dynamical effects produced by adding 
charcoal to fertile soils. This turns to be a nonlinearly coupled system of 
deterministic ODEs which behaves well mathematically. The main task is 
to understand to which extent charcoal (in its biochar form) is able to lock 
up carbon in soils. Our results are only preliminary in the sense that we are 
not solving here the CO2 sequestration problem. Rather we are setting up 
a modeling framework where this can be tackled by means of mathematical 
tools. 

The paper is organized as follows: In section|2]we describe mathematically 
chemical reactions in homogeneous media (here: fertile soils) and propose 
a first model based on differential equations. We prove in section [3] that 
our model is well-posed in the sense of Hadamard and perform a stability 
analysis of the interesting steady states. We illustrate the behavior of the 
profiles of the active concentrations and parameter effects in section |5j The 
effects observed regarding the addition of charcoal to soils are summarized in 
section Appendix A contains a discussion of the equilibria and stability 
of a reaction sub-block. 

We hope that our paper will bring the attention of the mathematical mod- 



eling community on the biochar issue - note that, cf. section 5.2, there are 



many open modeling aspects that would deserve a careful multi-disciplinary 
attention. 
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2. Modeling chemical reactions in homogeneous fertile soils 

2.1. What happens if charcoal is added to soil? 

In this section we provide a simple model for the chemical reactions tak- 
ing place in charcoal-enriched soil. We model only those processes that are 
relevant to carbon dioxide emission: the break down of soil organic matter 
and charcoal by microbes and the subsequent release of carbon dioxide, the 
reproduction and death of the microbes, and the effect of charcoal on soil 
fertility. 

We denote the species appearing in the chemical reactions by 
CO2 carbon dioxide, 

Ch charcoal (artificially added to the soil), 
Om soil organic matter (natural soil carbon), 
M microbes. 

Note that we do not distinguish between different types of soil organic mat- 
ter (litter, recalcitrant organic matter, humus, etc.). Also we only consider 
heterotrophic microbes, i.e., those that use organic carbon for growth. 

Microbes in the soil break down the organic matter and charcoal (this is 
called mineralization), releasing the carbon, which then combines with oxy- 
gen to form carbon dioxide. Experimental evidence indicates that generally 
there is no shortage of oxygen in the soil. Having this mind we assume that 
oxygen is present everywhere in equal amounts and thus it enters our model 
as a parameter. We model the complex system of mineralization processes 
by means of the following chemical reactions mechanism: 

Om ^ nC0 2 , (2.2) 
Ch % C0 2 , (2.3) 

where n > is taken as a constant. The reaction "constants" k\ and ki 
depend generally on the concentration of microbes, i.e, 

1* = h(M). 

Here we assume that, as functions, these reaction constants increase if the 
concentration of microbes increases. Note that, in general, the reaction con- 
stants can also depend on other effects (like the concentration of phenolics 

4 



(2.1) 



in the soil) , but for the sake of keeping things simple we do not include these 
in our model. 

The microbes need organic matter and oxygen to reproduce. Since we 
assumed that there is an abundance of oxygen, we can model the reproduction 
of microbes by means of 

M + 50m + 1)M, (2.4) 

where 5, fi > are constants. In general the reaction constant k% might de- 
pend on the fertility of the soil, which in turn depends on the amount of 
charcoal in the soil. For our theoretical investigations, we neglect the inter- 
mediate step and assume directly that ks depends on the amount of charcoal, 
^3 — ks(Ch), and that k 3 increases with charcoal concentration. However, 
note that the fertility of the soil contains so much in situ information that it 
cannot be neglected in the practical design of a CO2 sequestration scenario 
or if one wants to understand why terra preta (or 'black earth' ) is so fertile. 
Furthermore, in practice k$ depends on many other factors, e.g., tempera- 
ture, moisture, soil type, but we assume that these are all constant and so 
they do not appear explicitly in our model. 

We model the death of microbes by the chemical reaction 

M % rjOm, (2.5) 

where rj > is a constant. 

2.2. Basics of chemical kinetics 

We denote the concentration of species A at time t by [A] (t), e.g., [CO 2] (t) 
is the concentration of CO2 in the soil at time t. In order to derive evolution 
equations for the species concentrations we use the simple reaction ansatz, 
see, e.g., [6|. This assumption essentially states that if our set of reactions is 
given by the mechanism 

n n 
^2 a ij-Ai "A ^ PijAu J = 1, ••• ,777,, (2.6) 
i=l i=l 

where n e N denotes the number of species A i7 m e N denotes the number 
of chemical reactions, and otij,fiij G M + are stoichiometric coefficients, kj 



5 



reaction constants, then the elementary reaction rates are given by 

n 

Tj (A h A 2 ,..., A n ) := k 3 H[Ai] ai > . (2.7) 



Balancing the mass of the active species A4, we easily derive the evolution 
equations for the concentrations [At], viz. 



d 
dt 



[Ai] = y^ y (Pij - a ij) r J (A, A, • • • , An) 



i = l, 



n. 



(2.8) 



Before applying this methodology to (2.2)-(2.5), we introduce a new no- 
tation, see Table 2.1[ which is more convenient for the analysis. For the sake 
of readability and clarity, we use both notations throughout this paper. 



Ml 


[Om] 




[M] 




[Ch] 


U4 


[C0 2 ] 



Table 2.1: Alternative notation for the concentrations. 



Remark 2.9. (Restriction to spatially homogeneous soils) Within the frame- 
work of this paper, we consider a "continuously stirred tank reactor" case, a 
scenario intensively used in chemical engineering. See, e.g., J6]. In terms of 
soils, this means that we focus our modeling on a single space location, where 
the measurements are made, and we follow how the information "flows" over 
physically-important timescales. To this end, we assume the soil to be homo- 
geneous in the sense that no spatial substructures ( typically called microstruc- 
tures) appear, i.e., all soil components (gravel, sand, solid nutrients, water, 
etc) are well-mixed. We postpone for later the study of the more realistic case 
when the soil heterogeneities will be explicitly taken into account in terms of 
porosities, tortuosities, permeabilities very much in the spirit of JW ( general 
theory of flows in porous media), |2l El Ef (multiscale approaches to the 
chemical corrosion of concrete, smoldering combustion and plant growth, re- 
spectively), /IT) / (accumulation of cadmium in plants). Also, at a later stage 
it would be interesting to study the effect of the charcoal's platelet-like mi- 
crostructure (see Figure \w on the efficiency of adsorption and desorption 
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of the nutrients. Most likely this would lead to a two-scale ODE system in- 
timately coupled with evolution equations for the transport and storage of 
nutrients. 



Applying the simple reaction ansatz to (2.2)— (2.5), and assuming addi- 



tionally that the system has a constant source s > of organic matter, yields 
the nonlinear coupled system of ODEs 



d_ 

dt 
d_ 

dt 
d 

dt 
d 

dt 



u± = -ki(u 2 )ui - 5k 3 (u 3 )u 2 u{ + r]k A u 2 + s, 

u 2 = nkz{uz)u 2 u{ - k 4 u 2 , 

u 3 = -k 2 {u 2 )u 3 , 

m 4 = nki(u 2 )ui + k 2 (u 2 )u 3 . 



(2.10) 
(2.11) 
(2.12) 
(2.13) 



The source s can be thought of as organic matter entering the soil from the 
surface in the form of dead leaves, plants, etc. This system also requires 
initial conditions. Their role is to incorporate the type of soil. Throughout 



the rest of this paper we study the system (2.10)-(2.13) 



3. Mathematical analysis of the system (2.10)— (2.13) 



We start by making some assumptions on the model parameters entering 



(2.10)-(2.13). These assumptions will be of a technical nature and will be 



used to prove global existence of positive and bounded concentrations Ui and 
to study the steady states of this nonlinear ODE system. 



3.1. Restrictions on the model parameters 
We assume that 

5 > 1. 



(3.1) 



Assumption (3.1 ), together with the assumptions given below on the constitu- 



tive functions fcj, ensure that the right-hand side of the system (2.10)-(2.13) 



is Lipschitz continuous, which guarantees that our ODE system admits a 
unique local classical solution. 

In addition to assuming that S,r],fJL,n > 0, we assume that 



8 > Tjjl. 



(3.2) 
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This condition is used in Section 3.3 to show that the solution to (2.10)-(2.13) 
does not blow-up in finite time. 

Since the ki are reaction constants, we assume that they satisfy ki > for 
i G {1,2,3,4}. Note however that ki are nearly never true constants; they 
often incorporate a certain dependence on important physical/environmental 
quantities (here: spatial location, temperature, soil fertility, oxygen content, 
water content, etc). Here we take k& to be constant and assume that the 
functions ki : R — > (0, oo), j £ {1,2, 3}, are Lipschitz continuous and strictly 
increasing. For example, ki being strictly increasing means that an increase 
of microbes in the soil leads to an increase in the rate of break down of 
organic matter. 

Finally, we assume that the initial concentrations are positive and bounded, 
i.e. Ui(0) = u°i e [0, oo), % e {1, 2, 3, 4}. 

3. 2. Positivity of concentrations 

In this section we show that the concentrations ui,u 2 ,u 3 ,U4 are nonneg- 
ative for all times if their initial values are nonnegative. It suffices to show 
for each i e {1, 2, 3,4} that if Ui = and Uj > for all j ^ i, then iti > 0. 
This turns to be a trivial exercise: 

ui(0, u 2 ,uz, Ui) = r}kiU 2 + s > 0, 
u 2 (wi, 0,u 3 ,u 4 ) = 0, 
u 3 (mi,u 2 ,0,u 4 ) = 0, 

m 4 (mi, u 2) M3, 0) = nki(u 2 )ui + k 2 (u 2 )u 3 > 0. 

3.3. L°° bounds on concentrations 

We prove that the concentrations Ui do not blow-up in finite time. Fix 
arbitrary initial conditions u®. Then, based on the result of Section 3.2, we 
can assume that Ui > for all z = 1, 2, 3, 4. 



From the positivity of Ui and fcj, it follows immediately from ( 2.12[ ) that 



l M 3| 



(3.3) 



Adding equation (2.10) to 77 times equation (2.11) gives 
d 



dt 



(ui + 7]u 2 ) = —ki(u 2 )ui - (5 - r]fi)k 3 (u 3 )u 2 u 1 + s < s. (3.4) 
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The inequality (3.4) follows from (3.2) and the positivity of the ki and ttj. 
From (3.4) we conclude that u\ and U2 satisfy L°° bounds on any finite time 



interval. The numerics suggest that this bound is independent of the length 
of this time interval, but we do not need this here; see section |5j 

Relying on the L°° bounds on Ui for i 6 {1,2,3} on any finite time 



interval [0, r], we can bound the right-hand side of (2.13) by a constant 
C(t). Integration yields the bound 



u A {t) < C{T)t + u\ 



(3.5) 



for all t e [0, r], which immediately gives a bound on u± on any time interval 

[o,4 



3.4- Well-posedness 

Based on the positivity and the L c 
with the Lipschitz continuity of the right-hand side of (2.10)-(2.13), we recall 



bounds on concentrations, together 



classical ODE theory (see [121 [13], e.g.) to prove the following result: 

Theorem 3.6. (Global solvability). Assume that the assumptions stated in 
section 3.1 hold. Then for any set of initial conditions iti(0) = > 0, 
the system ( |2.1Q ) — ( 2.13 ) has a unique classical solution Ui : [0, oo) — >■ E, 
ie {1,2,3,4}. 

Furthermore, a Gronwall-like argument can be employed to show that 
this classical solution depends continuously on the initial data and all model 
parameters. Since this argument is rather standard, we omit to show it here. 



3.5. Equilibria and stability of the system (2.10)-(2.12) 



First note that M4 does not appear in the right hand side of (2.10 )-( 2.13 ). 



Hence equation (2.13) decouples from the system, in the sense that we do 



not need (2.13) to solve the subsystem (2.10 )-(2.12 ). Having this in mind, 



it is sufficient to study the equilibria of the reduced system (2.10 )-( 2.12 ). 



The reader is referred to Appendix A for a discussion of the equilibria and 



stability of the reaction block given by (2.4) and (2.5). For basic notions of 
dynamical system^} see [IE] , e.g. 



4 Dynamical systems theory proved to be very successful in a series of cases arising in 
biology and ecology; compare for instance [21 [15] and references cited therein. We expect 
therefore that dynamical systems delivers results in the case of biochar research as well. 
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We first search for the equilibria of the decoupled system given by (2.10), 



(2.11) and (2.12). By equating the right-hand side of (2.12) to zero, it follows 
that W3 = 0. By substituting this into equations (2.10) and (2.11) we obtain 








—ki(u 2 )ui - 5k^{Q)u 2 u[ + r]k 4 u 2 + s, 
(/x/c 3 (0)?4 - h)u 2 . 



(3.7) 
(3.8) 



For convenience we write k% instead of fc 3 (0) in the remainder of this section. 
Equation (3.8) is satisfied if and only if either 

u 2 = 0, or 



Ux 



/ik 3 



Let us treat the two cases separately: 



Case (3.9): It immediately follows from (3.7) that U\ = s/ki(0). 



Case (3.10): By inserting (3.10) in (3.7) we get 

k 4 

= -k 1 (u 2 )C 2 (5 - rjfi)u 2 + s. 



(3.9) 
(3.10) 



(3.11) 



The right-hand side of (3.11) is strictly decreasing as a function of u 2 . 
Hence it has at most one solution u 2 . A necessary condition for the 
existence of such a solution is that the right-hand side is nonnegative 
for u 2 = 0. This is the case when 

s > h(0)C 2 . (3.12) 

From now on we assume that the ki and the parameters 5, i], fi are 
chosen such that (3.11) has a solution whenever (3.12) holds. We will 
call this solution u\. For example, a solution exists if (3.12) holds and 

5 — TJ/jL > 0. 

Therefore, depending on the parameter s, we have one or two equilibrium 
points: If s < £4 (0)6*2, then we have only one equilibrium point u\ given by 



u\ := (ui,u 2 ,u 3 ) 



fci(0) 



,0,0 



(3.13) 
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If s > ki(0)C 2 , we have the additional equilibrium point u 2 given by 



U 2 e '■= (ui,U 2 ,U 3 ) 



(3.14) 



where u 2 satisfies (3.11). Therefore s = ki(0)C 2 is a bifurcation point. 

To test the stability of the equilibrium points u\ and u 2 e1 we linearize 



the system (2.10), (2.11), (2.12). Let J denote the Jacobian matrix of this 



system. A brief calculation shows that 



Mo) 
o 





kj(0) 



fci(O) 










Mo). 



(3.15) 



The eigenvalues of J(u\) are given by the entries on the diagonal. The eigen- 
values — fci(0) and — /c 2 (0) are negative, whereas the sign of the third eigen- 
value changes from negative to positive as s passes the bifurcation point. So 
u\ is asymptotically stable if s < —ki(0)C2 and is unstable if s > —ki(0)C2- 
We follow the same procedure for u 2 e . First we obtain 



M u %) ~S 2 k 3 u* 2 C d 2 
8nhu* 2 C s 2 - 1 




5-1 








Let us denote the 2x2 upper-left block of J{u 2 ) by 



6 fc 3(°) u „*' 
/i fc 3 (0) fi ' 4 "2 

fc 3 (0) /i4 "2 

-h(u* 2 ) 

(3." 



6) 



A 1 A 2 
A 3 



Note that Ai, A 2 < and A 3 > 0. Therefore the eigenvalues of Jiu 2 ) are 



-k 2 (u* 2 ), 



A, 



+ 



+ A 2 A 3l and 



A x 



+ A 2 A 3 . 



Since A\ < and A 2 A 3 < 0, the real parts of all the three eigenvalues are 
negative, which proves that u 2 is asymptotically stable. 

In summary, for each s > there is one stable equilibrium of the decou- 



pled system (2.10), (2.11), (2.12). Depending on the size of the source s, this 
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equilibrium is either given by (3.13) or by (3.14). Note that the full system 
(2.10 )— (2.13 ) does not have any equilibrium points since > (unless s — 0, 
in which case Ui = for alH = 1, 2, 3, 4 is an equilibrium). 



4. Nondimensionalisation 



Before solving the system numerically, we rescale it (very much in the 
spirit of [1TJ ) . We consider the following scalings for the time, concentrations, 
and reaction rates: t = rt, where r is the reference time, Ui = Uiiti, where 
Ui is the reference concentration of species i, and ki = Kiki, where Ki is the 



reference reaction constant. Substituting these into equations (2.10)-(2.13) 
gives 

-TifciWi - T 2 k 3 U 5 l U 2 + T 3 k A U 2 + T 4 S 



d _ 



—Mi 

dt 
d „ 

—Uo 

dt 



£ U3 



d 



dt 



~U4 



r b k\u[u 2 - r 6 k 4 u 2 , 
-r 7 k 2 u 3 , 

TshUx + Tgk 2 U 3 , 



(4.1) 



where r a ,a G {1,2, .. . ,9}, denote the characteristic time scales. Table 4.1 
lists their dependence on the reference constants. 



Characteristic 


Typical size 


time scale 




n 




T 2 


T5K 3 Ui-"U 2 


T 3 


Tr]K A U 2 Ur l 


TA 


tut 1 


T 5 


r^K 3 Ul 




tK 4 


T7 


tK 2 


T8 




T9 


tK 2 U 3 Ui L 



Table 4.1: List of characteristic time scales and their expected size. 



12 



parameter 


value 


ref. constant 


value 


unit 


a x 


1 


c 2 


1 


mol m ° 


a 2 


1 




0.01 


s 1 


a 3 


1.9 


T r 


10 


— 1 

s 


bi 


1 


K 3 


1 


m 36 mol _<5 s _1 


b 2 


1 




0.1 


s" 1 


b 3 


0.1 




1 


mol m -3 


W 


1 


u 2 


1 


mol m -3 


V 


10 




1 


mol m -3 


fx 


1 


u 4 


10 3 


mol m~ 3 


5 


10 


s 


0.02 


mol m -3 s _1 


n 


10 


T 


1 


s 



Table 5.1: Parameter values (first two columns) and reference constant values (last three 
columns) used for the simulations. 



5. Numerical simulation of the system (2.10)— (2.13) 



Here we illustrate numerically the behavior of the solution to our ODE 
system and test the effects of the various parameters. The main interest lies 
in predicting how the emission of C0 2 into the atmosphere changes if we put 
charcoal in the soil pQ. We make such predictions for different parameter 
values. 

We start by choosing the following linear constitutive functions for the 
reaction rates: 



k\{u 2 ) — a\u 2 + &i, k 2 (u 2 ) = a 2 u 2 + b 2 , 
k 3 (u3) = a 3 u 3 + b 3 , k 4 = b 4 . 



(5.1) 



The coefficient values are given in Table 5.1, along with the reference values 
C 2 , Ki, Ui, s and r, and the parameters rj,fj,,8,n. Note that we have chosen 



the parameters so that S = rjfx. From (3.10) it follows that 



C 2 



fiK 3 b 3 



We take s = 2if 1 a 1 C 2 so that u 2 e (see (3.14)) is a stable equilibrium point of 



the reduced system (2.10)-(2.12) 
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5.1. Searching for the effects of adding charcoal to soil 

First we simulate what happens to the equilibrium state u 2 e when charcoal 
is added to the soil. This corresponds to the following initial condition: 

(«i,« a ,«s,«*)U=^, Kiai c 2 U 2 ^°)- (5 - 2) 

We have exploited the special form of k\ to calculate u\ explicitly. 

Although the simulation is carried out for the dimensionless Ui, we will 
refer to them by [Om], [M], \Ch] and [CO2] for clarity. Figure [2] shows the 
simulation results. We see here various interesting phenomena: 

• The concentrations [Om] and [M] change on a short time scale (0(t) = 
0.1). Essentially this is because their time derivatives depend on [Ch] 
through k 3 . 

• On a long time scale (0(t) = 100), [Ch] decreases exponentially fast to 
0. Therefore [Om] and [M] converge back to their initial, equilibrium 
values. 

• On the same long time scale, [CO2] increases linearly. 




Figure 2: These figures show short-time (left) and long-time (right) behavior of the system 
initially at equilibrium. Charcoal is added at time t = 0. The graph of [CO2] is omitted 
in the left picture, because it would be too close to the t-axis to see anything interesting. 
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Now we turn our attention to CO2 emission. In Figure [3] we compare the 
case in which we put charcoal in the soil at time t — 0, which corresponds to 



initial condition (5.2), with the case in which the soil contains no charcoal, 



which corresponds to initial condition 



[Wl,W2,«3,W4j 



!i=0 



C2 s - K x C 2 b x 
£V K iai C 2 U 2 ,U,U 



(5.3) 



Initial condition (5.3) implies that [Om], [M] and [Ch] are constant in time 
and that [CO2] is linear. 




O 

£ 0.5 



—no [Ch] 
— with [Ch] 



300 2000 4000 6000 

t 



Figure 3: These figures show short-time (left) and long-time (right) emission of C02- We 
compare the case in which charcoal is added to the soil (in red) with the case when the 
soil does not contain charcoal (in blue). 

From Figure [3] we see that, on the short time scale Oil) = 100, putting 
charcoal in the soil reduces the CO2 emission. On the other hand, on a 
longer time scale {Oil) = 1000), we see that putting charcoal in the soil has 
no effect on the CO2 emission; the charcoal concentration [Ch] tends to zero, 
[Om] and [M] return to their equilibrium values, and the [CO2] concentration 
tends to the linear profile. 

A natural question that arises is whether we see similar effects if we 
increase the amount of charcoal that we put in the ground initially. We can 
simulate this by redoing the previous simulation, but now with U3 = 10 mol 
m -3 so that the amount of charcoal is ten times as much. The results are 
shown in Figure |4} On the long time scale the behaviour is similar to before. 
This is remarkable, because it means that the total amount of emitted CO2 
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hardly changes when ten times as much charcoal is put into the soil. On the 
short time scale we do see a difference: the rate of CO2 emission is slightly 
increased. 




Now we test the effect of some of the other parameters on the CO2 emis- 
sion. Figure [5] shows the results of doing the simulation with the values in 
Table 5.1 but now with K2 decreased by a factor of ten, i.e. K2 = 10~ 4 s _1 . 
This corresponds to a slower breakdown of the charcoal by the microbes. 
The qualitative behaviour is the same as Figure [3j but now the correspond- 
ing time scales are larger. This is because the charcoal is in the system for 
longer, and hence it takes longer to reach equilibrium. 

Figure [6] is the result of repeating the simulation shown in Figure [3] again, 
but now with the nonlinear constitutive function k\(u2) = aiu^ 00 + b±, which 
corresponds to a very fast breakdown of the organic matter by the microbes. 
The other parameters are the same as in Table 5.1 and equation (5.1 ). Note 
however that changing k\ changes the equilibrium u 2 e and so changes the 
initial condition. 

We observe that initially the rate of CO2 emission is much higher. This is 
because initially the number of microbes increases due to the presence of the 
charcoal, which speeds up the mineralization process tremendously because 
of the exponential constitutive function. On the long time scale, however, 
adding charcoal to the soil makes little difference. 
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Figure 6: The results of re-running the simulation shown in Figure [3] but with the 
nonlinear constitutive function fci(£s2) = aiu^ 00 + b%. 

5.2. Summary of our results and open problems 

For a rather large range of parameter values, our simulations clearly indi- 
cate that the short-time behaviour of our system can be significantly different 
from the long-time behaviour. Therefore, when testing experimentally the 
effect of adding charcoal to soil on CO2 emission, it could be important not 
to make judgements based solely on short-time data. 
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Our numerical results also suggest that the equilibrium v? e is globally 
asymptotically stable. 

In all of our simulations the long time CO2 emission seems to be inde- 
pendent of the addition of charcoal to the soil. This curious effect requires 
further study and possibly a refinement of our model. 

We are well-aware that we have a large set of parameters and reference 
constants that are difficult to relate to available experimental data. Therefore 
further simulations and comparison with experiments are required. This 
would naturally lead to a better control of the size of the characteristic time 
scales and potentially allow for improved predictions on C0 2 sequestration. 
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Appendix A. Equilibria and stability of the reaction block given 
by pi} and gjjb 



In this appendix we consider a subsystem of (2.10)-(2.13) that corre- 
sponds to reactions (2.4) and (2.5) (without the presence of charcoal, carbon 
dioxide or a source of organic matter). The reason for studying this sub- 
system is that it gives us a physical reason for imposing (3.2). Moreover, 
this subsystem turns out to dominate the short time behaviour of the whole 
system. 

Substituting s = 0, k x = and [Ch] = into fl2.10p and ( ]2.11| ) gives 
d 



dt 



[Om] = -5k 3 (0)[M}[Om} d + rik A [M], 



d 
dt 



(A.l) 



[M] = fik 3 (0)[M][Om} 5 - k A [M}. 



In the rest of this subsection we write k 3 instead of k 3 (0) for brevity. 

Figure A. 7 shows a sketch of the phase plane corresponding to (A.l). 
Note that 



d 
dt 



[Om] = [M] = or [Om] 



d_ 
dt 



[M] = & [M] = or [Om] 



5k 3 

k^_ 
fxk 3 



=:C 2 . 



(A.2) 



From (A.2) we see that ([Om],[M]) = (c, 0) is an equilibrium solution of 
flAlj ) for all ceR. If C\ — C 2 , then so is ([Om], [M]) = (Ci, c) for all ceR. 

To determine the stability of the first equilibria, ([Om], [M]) = (c, 0), we 
compute the Jacobian matrix corresponding to the system ( |A.l ): 



-5 2 [M][Om 
6jj[M}[Om 




- [Om] 



(A.3) 



From (A.3) it easily follows that the equilibria ([Om], [M]) = (c, 0) are stable 
if c < C 2 . 

Now we consider the boundedness of the trajectories. We consider three 
cases: C\ < C 2 , C\ = C 2 and C\ > C 2 (sketches of the corresponding phase 
planes are given in Fig. A.7). These cases correspond to: 
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Figure A. 7: Sketches of the phase plane corresponding to ( A.l ), depending on whether C\ 
is bigger or smaller than C2 (see (A.3) for their definitions). Recall that u\ — [Om] and 
u 2 = [M]. 



(5 > rjfi) : From the phase field analysis, we expect the solution of (A.l) to 
be bounded for all initial conditions. 



(5 = riii) : From ( A.2[ ), we see that we have more equilibrium points, which 
are given by [Om] = C\ = C2 and [M] e ffi. arbitrary. These equilibrium 
points are stably if and only if [M] > 0. 

(5 < riii) : From the phase field analysis, we expect the solution to blow up 
for most initial conditions. 



Therefore a sufficient condition for a solution of the reduced system (A.l) to 
be finite in time is 

5 > riii. (A.4) 



This is the same as our assumption (3.2) for the whole system. Equality 



in (A.4) would mean that the amount of organic matter that is converted 



into microbes by reaction (2.4) is equal to one over the amount of microbes 
that is converted into organic matter by reaction (2.5). This means that 



[Om] + ri[M] is conserved. Indeed, one sees immediately from (A.l) that 



d 
dl 



{[Om]+ri[M]) = 0. 



This quantity [Om] + rj[M] was also useful for proving L°° bounds for the 
whole system. See equation ( |3.4[ ). 
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